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EXECUTIVE  SUMMARY 


Maximum  Likelihood  Detection  Algorithm 

The  description  of  a  maximum  likelihood  algorithm  to  detect  moving  targets  in  electro-optic 
data  is  presented.  The  algorithm  is  evaluated  in  terms  of  the  probabilities  of  false  alarm  and  detec¬ 
tion.  A  comparison  of  theoretical  and  experimental  probability  distributions  for  single  normalized 
pixels  shows  good  agreement.  Similarly,  a  comparison  of  theoretical  and  experimented  false  alarm 
probabilities  also  shows  good  agreement.  These  results  validate  using  theoretical  models  to  predict 
algorithm  performance.  Several  detection  examples  are  shown. 

The  data  sets  compared  against  theory  are  obtained  from  different  sensor  types,  wavebands, 
and  clutter  backgrounds.  The  Experimental  Test  System  sensor  provides  visible  band  data  from 
a  staring  sensor.  The  background  clutter  in  these  data  sets  consists  of  stars  and  other  deep  space 
objects.  The  Infrared  Measurements  Sensor  (IRMS)  provides  long  and  medium  wavelength  infrared 
data  from  a  scanning  sensor.  Background  clutter  in  the  IRMS  data  sets  consists  of  sky,  mountains, 
hills,  and  desert. 

Binary  Integration  Algorithm  and  Architecture 

A  binary  integration  version  of  this  algorithm  is  described  and  evaluated  in  terms  of  false 
alarm  and  detection  probabilities.  This  version  is  suboptimum  and  is  compared  with  the  optimum 
algorithm  to  determine  the  performance  loss.  A  processing  architecture  concept  is  also  described. 

Conclusion 

The  model  used  in  the  maximum  likelihood  algorithm  development  is  shown  to  closely  agree 
with  experimental  data  when  the  clutter  background  has  temporally  stationary  statistics.  An 
accurate  model  allows  the  algorithm  performance  (probabilities  of  detection  and  false  alarm)  to 
be  precisely  predicted.  Potential  areas  for  future  work  include:  (1)  development  of  processing 
architectures  and  efficient  algorithm  implementations  for  fast  computation,  (2)  inclusion  of  target 
signal  (signal-to-noise  ratio)  statistics  in  the  model  to  investigate  improved  algorithms,  and  (3) 
investigation  of  frame  registration  issues  and  techniques. 

Maximum  Likelihood  Algorithm  Derivation 

An  algorithm  to  detect  moving  targets  in  electro-optic  data  is  derived  for  the  case  where  the 
data  samples  have  Gaussian  (normal)  distributions  that  are  temporally  stationary  (frame  to  frame) 
and  spatially  nonstationary  (pixel  to  pixel).  In  other  words,  the  means  and  variances  describing 
the  data  samples  are  assumed  to  vary  among  pixels  according  to  the  background  scene  but  remain 
constant  over  time.  Additionally,  these  underlying  means  and  variances  are  presumed  unknown. 
Targets  are  assumed  to  have  an  arbitrary  distribution. 


in 


Single  Pixel  Probability  Distributions 

The  probability  distribution  of  a  single  normalized  pixel  (used  in  the  likelihood  function)  is 
derived.  This  distribution  is  used  elsewhere  to  derive  the  theoretical  probabilities  of  false  alarm 
and  detection. 

In  the  case  of  noise-only  (no  target),  the  probability  distribution  is  shown  to  be  an  instance  of 
a  beta  distribution.  With  a  target  present,  the  probability  distribution  is  a  noncentral  beta  distri¬ 
bution.  These  distributions  depend  on  the  background  clutter  scene  only  through  the  noncentrality 
parameter,  which  is  proportional  to  the  signal-to-noise  ratio. 

Likelihood  Function  Statistics 

The  probabilities  of  false  alarm  and  detection  are  derived  for  the  maximum  likelihood  detec¬ 
tion  algorithm.  Evaluation  of  these  probabilities  is  accomplished  through  numerical  integration. 
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1.  MAXIMUM  LIKELIHOOD  DETECTION  ALGORITHM 


An  algorithm  to  detect  moving  targets  in  electro-optic  data  is  developed  and  analyzed.  In 
such  data  moving  targets  are  distinguished  by  their  motion  relative  to  a  fixed  background  scene 
(clutter).  Unresolved  targets  may  appear  as  a  “streak”  in  a  single  image  frame  collected  with  a 
staring  sensor  or  as  a  dot  in  a  single  image  frame  from  a  scanning  sensor.  A  sequence  of  images 
can  be  used  to  distinguish  between  clutter  and  targets. 

A  maximum  likelihood  algorithm  is  derived  to  detect  a  target  with  an  arbitrary  signal  distri¬ 
bution  against  unknown  clutter.  A  consequence  of  this  algorithm  is  that  detection  thresholds  may 
be  set  according  to  a  given  false  alarm  rate  without  knowledge  of  the  background  clutter  scene. 
Experimental  false  alarm  performance  is  compared  against  the  theoretical  predictions  and  found 
to  closely  agree. 

The  data  used  to  test  this  algorithm  are  obtained  from  two  sensors  developed  at  Lincoln 
Laboratory.  The  Experimental  Test  System  (ETS)  sensor  is  a  visible  band  staring  sensor  that  is 
ground-based  and  typically  used  to  detect  moving  objects  against  a  nighttime  stellar  background. 
In  these  data  clutter  consists  of  stars  and  other  deep  space  objects  as  well  as  atmospheric  effects, 
notably  clouds.  Figure  1  shows  an  unprocessed  image  from  the  visible  band  test  data  sequence.  A 
pseudocolor  mapping  indicates  intensity. 

The  Infrared  Measurements  Sensor  (IRMS)  is  a  scanning  sensor  that  collects  both  long  wave¬ 
length  infrared  (LWIR,  8  to  12  /im)  and  medium  wavelength  infrared  (MWIR,  3  to  5  /im)  data. 
This  sensor  is  typically  used  with  a  clutter  background  of  sky,  mountains,  desert,  and  forest.  Both 
sensor  types  work  well  with  this  algorithm.  Figures  2  and  3  show  unprocessed  images  from  three 
of  the  IR  test  data  sequences.  Again,  a  pseudocolor  mapping  indicates  intensity. 

The  technical  discussion  that  follows  is  divided  into  four  parts.  A  statement  of  the  maximum 
likelihood  detection  algorithm  is  given  first.  The  theoretical  probability  distribution  of  normalized 
data  is  second  with  a  comparison  of  the  theoretical  and  experimental  results.  The  theoretical 
probability  of  false  alarm  and  false  alarm  rate  are  third,  also  with  a  comparison  of  the  theoretical  and 
experimental  results.  The  final  section  presents  the  theoretical  detection  probability  and  includes 
several  detection  examples. 

1.1  Signal  Model  and  Algorithm 

The  maximum  likelihood  moving  target  detection  algorithm  is  derived  under  the  assump¬ 
tion  that  data  samples  collected  from  a  sequence  of  image  frames  are  independent,  normally  (i.e., 
Gaussian)  distributed  random  variables.  In  the  absence  of  targets,  these  samples  are  modeled  as 
stationary  in  time  but  nonstationary  in  space.  In  other  words,  the  means  and  variances  of  the  data 
samples  are  assumed  to  vary  from  pixel  to  pixel,  but  for  a  single  pixel  they  are  constant  over  time. 

The  motivation  for  an  independent  normal  statistics  model  is  based  on  assumptions  of  (1) 
independent  thermal  and  readout  sensor  noise  and  (2)  independent  arrival  of  photons  to  each 
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Figure  1.  Single  unprocessed  frame  of  ETS  Magellan  data 
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Figure  2.  Single  unprocessed  frame  of  IRMS  Longjump  data 
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Figure  3.  Single  unprocessed  frame  of  IRMS  Dugway  data 
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detector  (a  Poisson  random  variable  that  may  be  approximated  as  normal).  Assuming  stationarity 
over  time  implies  assuming  that  all  frames  in  a  sequence  are  spatially  registered  and  that  there  is 
no  moving  clutter.  Nonstationary  spatial  statistics  correspond  to  the  nonuniform  intensity  of  the 
background  scene. 

An  equivalent  viewpoint  is  to  consider  the  mean  for  a  given  pixel  to  be  the  image  intensity  in 
that  pixel  and  the  variance  to  correspond  to  the  (photon  plus  thermal)  noise  level.  The  indepen¬ 
dence  of  samples  then  corresponds  to  each  sample  having  an  independent  noise  component,  while 
the  means  correspond  to  the  imaged  scene. 

The  log-likelihood  function  L  (which  results  from  this  data  model)  follows,  where  rtJt*  is  the 
data  sample  in  row  i,  column  j,  frame  A;,  and  N  is  the  number  of  frames  available  (see  Appendix  A). 


£  =  £  d{Jtk  (1) 


(ij,k)es 

(riJ,k  fit j)2 

a-?. 

(2) 

= 

l  N- 1 

Tv 

n=0 

(3) 

ah  = 

1  N~' 

■*  ' _ n 

(4) 

The  exact  underlying  means  and  variances  are  /itj  and  cr^  and  are  presumed  unknown.  Target 
detection  is  accomplished  by  comparing  £  to  a  threshold. 

Set  S  is  the  set  of  space  (i,j)  and  time  ( k )  indices  that  corresponds  to  the  target  motion 
hypothesis.  As  an  example,  consider  a  target  that  moves  with  speed  v  pixels/s  in  the  (i,j)  plane 
at  an  angle  0 .  For  a  frame-to-frame  period  of  T  s,  set  5  is 


S  =  {(i  +  vTkcos(6),  j  +  vTk  sin(0),  k);  k  =  0, . . . ,  AT  —  1}  ,  (5) 

where  the  noninteger  coordinates  are  rounded  to  the  nearest  integer  values. 

1.2  Single  Pixel  Statistics 

The  term  dtj^  in  the  formula  for  the  log-likelihood  function  is  referred  to  as  the  “normalized” 
data.  Selection  of  detection  thresholds  and  the  resulting  probabilities  of  false  alarm  and  detection 
depend  on  the  statistics  of  d{jfk  and  C.  The  theoretical  distribution  jFj(x|At)  of  d{jk  has  been 
derived  and  shown  to  be  an  instance  of  the  beta  distribution  (see  Appendix  B). 
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=  Fb 


X 


(6) 


TV  -  1 


1  N-2 
2’  2 


It  is  interesting  to  note  that  the  statistics  of  d{jtk  depend  only  on  the  number  of  frames  processed 
and  not  on  the  terms  /it-j  and  related  to  the  specifics  of  the  background  clutter.  Thus  the 
effects  of  clutter  are  removed  from  determination  of  false  alarm  rates. 

Clutter  does  have  an  effect  on  the  single  pixel  statistics  when  targets  are  present.  The  pixel 
variance  cr2j  is  the  noise  component  of  the  signal-to-noise  ratio  (SNR)  and  is  related  to  the  pixel 
mean  /!{ j .  (Specifically,  if  the  data  values  were  the  number  of  photons  received,  then  the  mean  and 
variance  of  a  pixel  would  be  equal.)  It  is  shown  in  Appendix  B  that  an  instantaneous  signal-to-noise 
ratio  SNRj  =  s1  jo2  can  be  defined,  where  s  is  the  target  amplitude  component  of  the  unnormalized 
pixel  data,  and  a2  is  the  noise  power  component.  Then  the  conditional  distribution  for  a  single 
normalized  pixel  is  the  noncentral  beta  distribution 


Fd(x\N,SKRi) 


x 

N  -  1 


1  TV  -2 
2’  2 


A7  -  1 
TV 


(7) 


Figure  4  shows  the  theoretical  normalized  pixel  probability  distribution  for  various  values  of  SNRj 
[expressed  in  decibels,  lOlog(SNRj)]. 

A  comparison  of  theory  with  empirically  derived  statistics  can  be  used  to  validate  the  data 
model  and  the  use  of  the  resulting  maximum  likelihood  detection  algorithm.  Figures  5  through 
7  compare  the  theoretical  distribution  for  a  single  normalized  pixel  and  experimentally  derived 
distributions  for  a  variety  of  image  sequences.  These  comparisons  are  done  for  the  noise-only  case 
(no  target). 

Figure  5  compares  theory  and  experiment  for  15  frames  of  visible  band  data  obtained  from 
the  ETS  sensor,  which  is  a  ground-based  staring  sensor,  looking  at  a  background  of  stars.  Although 
the  stellar  background  is  spatially  nonstationary,  there  is  good  agreement  between  the  two  curves. 

Figure  6  compares  theory  and  experiment  for  10  frames  of  the  Longjump  LWIR  data  obtained 
from  the  IRMS  scanning  sensor,  which  is  looking  at  a  background  of  sky,  mountains,  and  hills. 
Again,  there  is  good  agreement. 

Figure  7  compares  theory  with  10  frames  of  LWIR  data  from  Dugway,  obtained  from  the 
IRMS  sensor.  Background  clutter  includes  sky,  mountains,  and  desert.  Departure  from  theory  is 
due  in  part  to  frame  registration  errors  but  may  include  other  transient  (time  varying)  interference 
as  well.  In  general,  data  that  are  not  temporally  stationary  will  disagree  with  this  theory. 

1.3  False  Alarm  Statistics 

Once  the  single  pixel  statistics  are  known,  it  is  possible  to  determine  the  theoretical  probability 
of  false  alarm  and  compare  this  theory  to  experimental  data.  The  probability  of  false  alarm  is 
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Figure  4-  Theoretical  normalized  pixel  probability  distribution  for  different  SNRj. 


Figure  5.  Single  normalized  pixel  probability  distribution  of  ETS  Magellan  data. 
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Figure  6.  Single  normalized  pixel  probability  distribution  of  IRMS  Longjump  data. 


Figure  7.  Single  normalized  pixel  probability  distribution  of  IRMS  Dugway  data . 
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defined  as  the  probability  that  the  log-likelihood  function  will  exceed  the  detection  threshold  when 
no  target  is  present.  System  false  alarm  rate  is  determined  by  combining  the  probability  of  false 
alarm  with  the  number  of  pixels,  the  number  of  velocity  hypotheses  tested,  and  the  data  collection 
time. 

Using  the  definitions 

Np  =  number  of  pixels  in  one  frame 
Nv  =  number  of  velocity  hypotheses 
N  =  number  of  image  frames  processed 
T  =  time  to  collect  one  image  frame 
t  =  detection  threshold 
R  =  false  alarm  rate, 

the  false  alarm  rate  may  be  determined  as  the  number  of  false  alarms  per  second  from  the  formula 


R  = 


NpNv?T{£>t\N} 

NT 


(8) 


The  probability  term  Pr{£  >  t\N}  is  computed  as  the  iV-fold  convolution  of  single  pixel  density 
functions  with  the  complement  of  a  unit  step  function  s(i). 


Pr{£>  t\N)  =  /<,*•••  * /«,*(!  -«)(<) 


(9) 


This  computation  is  implemented  using  numerical  integration  on  the  beta  density  functions  (see 
Appendix  C). 

As  an  example  of  the  false  alarm  rate,  testing  1,000  velocity  hypotheses  on  10  frames  of  IRMS 
data  results  in 


R 


(400  x  3480) x  1000  x  10"8 
10  x  1 


(10) 


=  1.392  false  alarms/s. 


(11) 


The  detection  threshold  t  is  assumed  to  be  chosen  so  that  the  probability  of  a  false  alarm  from  a 
single  target  position  and  velocity  test  is  10-8.  The  full  frame  of  the  IRMS  data  is  used  with  a 
corresponding  collection  rate  of  1  frame/s. 
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The  validity  of  the  false  alarm  rate  calculation  depends  on  the  accuracy  of  the  probability  of 
false  alarm  term  in  the  expression  for  R.  For  the  probability  of  false  alarm  term  to  show  agreement 
between  theory  and  experiment  at  a  level  of  10-8,  the  theoretical  and  experimental  probability 
distributions  (i.e.,  Pr{£  <  t\N}  =  1  -  Pr{£  >  t\N})  of  the  log-likelihood  function  should  agree  to 
eight  significant  digits. 

Theoretical  and  experimental  probabilities  of  false  alarm  are  compared  in  Figures  8  through 
10.  Experimental  probabilities  are  determined  by  testing  104  velocity  hypotheses  on  1.8  X  105  pixels 
of  the  visible  band  data  and  1.4  X  105  pixels  of  IR  data.  Thus  the  total  number  of  hypotheses  tested 
for  each  experiment  is  approximately  109.  For  each  data  set,  a  histogram  is  formed  from  the  values 
of  the  log-likelihood  function  for  these  trials.  The  probability  curve  is  then  easily  computed  from 
the  histogram. 

Figure  8  shows  false  alarm  statistics  for  15-frame  processing  of  the  ETS  Magellan  visible  band 
data  set,  and  Figure  9  shows  10-frame  processing  of  IRMS  Longjump  LWIR  data.  These  data  sets 
show  good  agreement  between  theory  and  experiment  (the  two  curves  are  difficult  to  distinguish!). 
Specifically,  the  two  curves  agree  at  probabilities  down  to  at  least  10~8.  Figure  10  shows  10-frame 
processing  of  IRMS  Dugway  LWIR  data. 


THRESHOLD /N 


Figure  8.  False  alarm  probability  for  ETS  Magellan  data . 
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PROBABILITY  OF  FALSE  ALARM  PROBABILITY  OF  FALSE  ALARM 
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Figure  9 .  False  alarm  probability  for  IRMS  Longjump  data . 
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Figure  10.  False  alarm  probability  for  IRMS  Dugway  data. 
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1.4  Detection  Probability  and  Examples 


As  noted  in  Appendix  C,  the  probability  of  detection  depends  on  the  SNR  of  the  normalized 
pixels  summed  in  the  velocity  hypothesis.  The  exact  SNR  distribution  depends  on  the  background 
clutter  and  is  difficult  to  specify. 

A  simple  case  of  the  SNR  distribution  for  uniform  background  clutter  is  discussed  in  Ap¬ 
pendix  C.  Assuming  that  the  dominant  noise  component  is  photon  noise,  the  noise  variance  (power) 
of  a  pixel  is  proportional  to  the  mean,  which  in  turn  corresponds  to  the  pixel  intensity.  Thus  the 
noise  component  of  the  SNR  depends  on  the  background  clutter  scene.  The  signal  (target)  compo¬ 
nent  is  also  derived  from  a  photon  count  and  in  the  absence  of  other  fluctuations  gives  a  noncentral 
chi-square  distribution.  For  uniform  intensity  background  clutter  the  SNR  distribution  can  be 
modeled  as  a  constant  times  a  chi-square  random  variable.  For  this  simple  case  the  detection  prob¬ 
abilities  depend  on  both  the  SNR  and  the  target  signal  power.  Other  SNR  models  also  depend  on 
the  clutter  intensity  distribution. 

As  an  example,  Figure  11  shows  the  theoretical  probability  of  detection  when  the  same  number 
of  target  photons  are  received  in  each  pixel,  and  a7  is  identical  across  all  pixels.  In  other  words, 
SNRj  is  assumed  constant.  In  this  case,  N  =  10  frames  are  assumed  to  be  processed. 

In  the  detection  examples  that  follow,  an  effective  signal-to-noise  ratio,  SNReff,  is  estimated 
from  the  data.  This  value  is  defined  as  the  value  of  SNRj,  which  has  a  normalized  pixel  mean 
(expected  value)  equal  to  the  sample  average  of  normalized  pixels  along  the  hypothesized  target 
path.  These  two  pixel  quantities  are  not  the  same,  because  the  expected  value  is  determined  for 
a  single  value  of  SNRj  =  s2/<r2,  and  the  sample  average  is  determined  from  data  where  both  s7 
and  a7  vary  among  pixels.  Thus  the  term  “effective  SNR”  is  used.  Figure  12  shows  the  normalized 
pixel  mean  as  a  function  of  SNRj  for  10-  and  15-frame  processing. 

Figure  13  shows  a  detection  example  for  the  ETS  Magellan  visible  band  data  set.  Approxi¬ 
mately  600  velocity  hypotheses  are  applied  to  these  data  (418  X  420  pixels/frame).  The  false  alarm 
rate  for  a  single  hypothesis  test  is  set  to  10~8,  and  the  expected  number  of  false  alarms  is  0.9. 
(The  target  motion  hypothesis  is  required  to  fit  entirely  within  the  data  in  this  experiment,  so  that 
not  all  combinations  of  position  and  velocity  are  used.)  The  estimated  values  of  SNRepp  for  these 
targets  are  10.3  and  15.3  dB. 

Figure  14  shows  a  detection  example  for  the  IRMS  Longjump  LWIR  data  set.  Approximately 
1,000  velocity  hypotheses  are  applied  to  a  400  X  400  subimage  of  the  data.  The  false  alarm  rate  for 
a  single  hypothesis  test  is  set  to  10-8,  and  the  expected  number  of  false  alarms  is  1.4.  (As  above, 
not  all  combinations  of  position  and  velocity  are  used.)  The  estimated  value  of  SNRe^*  for  this 
target  is  17.3  dB. 

Figure  15  shows  a  detection  example  for  the  IRMS  Dugway  LWIR  data  set.  Approximately 
1,800  velocity  hypotheses  are  applied  to  a  400  X  500  subimage  of  the  data.  The  false  alarm  rate  for 
a  single  hypothesis  test  is  set  to  10”8,  and  the  expected  number  of  false  alarms  is  1.4.  (Again,  not 
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Figure  11.  Theoretical  probability  of  detection  for  constant  SNRt. 
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Figure  12.  Expected  value  (mean)  of  normalized  pixels  as  a  function  of  SNR,. 
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all  combinations  of  position  and  velocity  are  used.)  The  upper  detection  is  a  target,  and  the  lower 
detection  is  a  false  alarm.  The  estimated  value  of  SNRe^  for  the  actual  target  is  9.5  dB. 

Target-related  false  alarms  can  occur  when  a  hypothesized  target  path  intersects  the  actual 
path  of  a  strong  target.  From  Equation  (1)  it  is  clear  that  a  single,  strong,  normalized  pixel  can 
cause  the  log-likelihood  function  C  to  exceed  a  detection  threshold.  One  way  to  reduce  the  false 
alarms  related  to  strong  targets  and  clutter  discretes  is  to  impose  a  model  of  the  target  statistics  in 
the  likelihood  function  derivation.  This  approach  has  the  disadvantage  of  depending  on  the  target 
and  clutter  models,  as  well  as  resulting  in  a  much  more  complex  formula  for  the  log-likelihood 
function.  A  different  approach  is  to  require  that  any  detection  also  has  at  least  a  certain  number  of 
strong  pixels.  This  concept  is  the  basis  for  the  binary  detection  algorithm  discussed  in  Section  2. 
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(a) 


Figure  IS.  Target  detection  for  ETS  Magellan  data,  (a)  Detected  target  locations,  (b) 
Extrapolated  target  motion. 
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(a) 


Figure  lj .  Target  detection  for  IRMS  Longjump  data,  (a)  Detected  target  locations,  (b) 

Extrapolated  target  motion. 
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(a) 


(b) 

Figure  15.  Target  detection  for  IRMS  Dugway  data,  (a)  Detected  target  locations,  (b) 
Extrapolated  target  motion. 


19 


2.  BINARY  INTEGRATION  ALGORITHM  AND  ARCHITECTURE 


For  many  applications,  the  algorithm  discussed  in  Section  1  is  computationally  intensive. 
This  section  describes  an  algorithm  for  a  binary  integration  approach  to  velocity  filtering.  Binary 
integration  is  of  interest  as  a  means  of  reducing  both  memory  and  computation  rate  requirements 
at  a  small  loss  (typically  less  than  2  dB)  in  processing  gain.  Also  presented  is  an  implementation 
architecture  that  is  suitable  for  a  small  applications  specific  integrated  circuit  (ASIC),  which  could 
be  attached  to  a  general  purpose  signal  processor. 

An  overview  and  performance  measures  of  binary  integration  are  given  in  Section  2.1.  The 
algorithm  is  described  mathematically  in  Section  2.2.  An  architecture  for  efficient  implementation 
is  described  in  Section  2.3.  Extension  of  this  algorithm  to  nonscanned  sensors,  where  targets  form 
streaks  within  an  image  frame,  is  discussed  in  Section  2.4. 

2.1  Algorithm  and  Performance 

Binary  integration,  a  technique  that  is  well  known  in  the  radar  community,  is  often  referred  to 
as  an  “M-out-of-N  detector”  [1,  2].  A  variation  of  binary  integration  is  also  used  in  the  Maximum 
Value  Projection  algorithm  implemented  on  the  Space-Based  Visible  processor  [3]. 

The  relationship  between  binary  integration  and  other  electro-optic  signal  processing  functions 
is  shown  in  Figure  16.  Typically,  sensor  data  may  be  corrected  (e.g.,  for  gamma  circumvention  in 
space  applications)  and,  in  the  case  of  scanning  sensors,  time-delay  integrated  to  improve  target 
SNR.  Most  moving  target  detection  algorithms  also  require  some  form  of  frame  registration.  The 
remainder,  data  normalization/quantization  and  velocity  filtering,  form  the  moving  target  detection 
algorithm.  Binary  integration  addresses  these  last  two  parts. 
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Figure  16.  Signal  processing  functions. 


The  binary  integration  approach  may  be  divided  into  four  steps: 

1.  Data  normalization 

2.  Binary  quantization 
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3.  Integration  (summation)  over  the  velocity  hypotheses 

4.  Threshold  detection  of  the  integration  results. 

Data  may  be  normalized  using  a  variety  of  techniques.  The  maximum  likelihood  algorithm  in 
Section  1  suggests  that  normalized  data  be  obtained  from  unprocessed  data  rIJt*  according 
to 


itjtk  — 

a? 

P'iJ  ~ 

1  N-l 

N  LI 

k=0 

ii 

<b 

1  yv_1 

N  ~  £'.j) 

A:=0 

(12) 


(13) 


(14) 


where  i,j  are  the  pixel  indices  and  k  is  the  image  frame  index.  Binary  quantized  data  6tJ  ^  are 
then  obtained  by  comparing  d{j^  with  a  fixed  threshold  t. 


0,  diJtk  <  t 

1  y  dtj^k  ^  ^ 

The  integration  step  is  also  called  velocity  filtering.  Samples  in  (*,j,  k)  are  integrated  along 
a  hypothesized  target  velocity  vector.  This  velocity  hypothesis  corresponds  to  the  velocity  filter. 

The  number  of  samples  used  in  integration  depends  on  sensor  type.  For  a  scanning  sensor,  a 
target  moves  less  than  one  pixel  during  the  time  the  target  is  scanned;  however,  the  total  scanning 
time  of  a  frame  is  typically  large  enough  so  that  a  target  moves  several  pixels  from  one  frame  to 
the  next.  Thus  a  target  streak  in  (z,j,  k)  appears  as  a  dotted  instead  of  solid  line.  In  this  case,  the 
integration  step  requires  summing  only  1  pixel/frame  for  spatially  unresolved  targets. 

Let  Sij  be  the  result  of  summing  (integrating)  over  the  data  for  a  target  that  has  initial 
position  z,j  and  frame-to-frame  velocity  components  u,  v  in  the  x,y  directions.  The  result  stj  is 
compared  with  a  threshold  in  order  to  achieve  a  detection.  (This  summation  and  comparison  is 
the  M-out-of-N  detector.)  The  quantity  S{j  is  computed  according  to  the  formula1 


aIn  general,  velocity  components  it,  v  are  floating  point  numbers;  indices  i  +  uk  and  j  +  vk  are  then 
taken  to  be  the  nearest  integer  values.  This  method  is  called  nearest  neighbor  interpolation. 
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(16) 


N- 1 

si,j  =  E  b<+uk,j+vk,k 

k=  0 


The  performance  of  the  binary  integration  algorithm  can  be  given  in  terms  of  the  probabilities 
of  detection  and  false  alarm  by  computing  the  probability  that  stiJ  meets  or  exceeds  the  threshold 
M  for  a  given  SNR.  Let  p  be  the  probability  that  a  single  normalized  pixel  exceeds  threshold  t, 

p  =  1  —  F<i(x|iV,  SNRj)  (17) 


=  1  - 


x 

N  -  1 


1  N- 2 

2’  2 


N  -  1 
N 


SNRj 


) 


(18) 


Using  the  constant  SNRj  example  discussed  in  Section  1,  the  probability  of  detection  Pp  is  obtained 
as  the  probability  of  a  binomial  random  variable  with  parameter  p  having  at  least  M  occurrences 
out  of  TV  trials. 


PD  =  Pr{s,j  >  M} 


(19) 


=  E  (  N)pmd-P)N- 

m—M  \  m  ) 


(20) 


The  corresponding  probability  of  false  alarm  Pfa  is  obtained  by  setting  SNRj  =  0. 

The  two  thresholds,  t  and  M ,  are  usually  chosen  to  maximize  the  Pp  for  a  particular  operating 
point  in  terms  of  SNRj  and  Pfa .  This  maximization  can  be  performed  by  computing,  for  each 
possible  value  of  M  =  0, . . . ,  TV,  the  value  of  t  that  gives  the  desired  Pfa •  The  (*,  M )  pair  to  be 
used  is  that  which  maximizes  the  value  of  Pp  for  the  required  SNRj. 

As  an  example,  consider  TV  =  10  frames  and  let  the  operating  point  for  this  optimization  be 
Pfa  =  10“8  and  SNRj  =  10.  The  maximization  of  Pp  yields  M  =  7  and  t  =  3.310.  Probabilities 
of  detection  and  false  alarm  can  be  plotted  for  a  constant  M  or  f. 

Figure  17  shows  the  probability  of  detection  as  a  function  of  <,  the  single  pixel  threshold,  for 
M  =  7  and  several  different  SNRj.  Figure  18  shows  the  false  alarm  probability. 

Figure  19  shows  the  probability  of  detection  as  a  function  of  Af,  for  t  =  3.310  and  several 
different  SNRj.  Figure  20  shows  the  false  alarm  probability. 

A  comparison  of  Figure  17  with  Figure  11  indicates  the  loss  in  SNR  for  the  binary  integration 
algorithm.  For  the  specified  Pfa  =  10”8,  the  binary  algorithm  requires  less  than  11-dB  SNRj  target 
to  achieve  the  same  detection  probability  as  the  full  precision  maximum  likelihood  algorithm  with 
a  10-dB  SNRj  target,  indicating  a  loss  less  than  1-dB.  This  loss  is  similar  to  that  obtained  when 
binary  integration  is  used  in  radar  processing  [1,  2]. 
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PROBABILITY  OF  FALSE  ALARM  <S  PROBABILITY  OF  DETECTION 
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are  17.  Probability  of  detection  for  binary  integrator ,  N  =  10  frames,  M  —  1. 


t 


Figure  18.  Probability  of  false  alarm  for  binary  integrator,  N  —  10  frames,  M  —  7. 
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Figure  19.  Probability  of  detection  for  binary  integrator  N  =  10  frames,  t  =  3.310. 


Figure  20.  Probability  of  false  alarm  for  binary  integrator,  N  —  10  frames,  t  —  3.310. 
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2.2  Efficient  Integration  Algorithm 

Computationally,  the  most  intensive  step  is  integration.  Efficiency  is  attained  by  reducing 
the  number  of  data  accesses  required  in  testing  a  given  velocity  hypothesis  on  many  possible  target 
positions. 

Let  Mx,  My  be  the  frame  size  in  the  x,  y  dimensions,  respectively.  There  are  ~  MxMy  possible 
initial  target  positions.  For  each  velocity  hypothesis,  integration  and  detection  of  N  data  frames 
can  be  accomplished  straightforwardly  by  summing  N  samples  (one  from  each  frame)  for  each 
possible  target  position.  This  method  requires  MxMyN  additions  or  comparisons  and  MxMyN 
data  accesses.  The  objective  of  the  algorithm  presented  here  is  to  reduce  the  number  of  data 
accesses  to  MxMy  and  to  perform  N  additions  with  each.  Thus  it  will  be  shown  that  the  MxMy 
possible  initial  target  positions  can  be  tested  at  an  average  cost  of  only  one  data  access  per  test.  The 
algorithm  derivation  that  follows  is  rather  mathematical  and  may  be  skipped  in  favor  of  reading 
the  step-by-step  description  given  at  the  end  of  this  section. 

The  first  step  is  to  pack  N  frames  of  binary  data  into  iV-bit  integers,  at|J, 


N-i 

a*\i  =  53 

k=0 


(21) 


where 

0  <  i  <  Mx  (22) 

0  <  j  <  My  .  (23) 

The  above  formula  for  can  be  written  as 
N—i 

=  5  v  7  (24) 

k=  0 


where  Bk(z)  is  bit  k  of  2. 

Now  consider  the  computation  of  Si+UJ+v. 


N—l 

i+u+uk,j+v+vk) 

k=0 


N 

y  Bk-i ( ai+uk,j+vk ) 
A:=l 


(25) 

(26) 
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It  is  important  to  note  that  N  —  1  of  the  data  accesses  required  (that  is,  indices  of  a,j)  are 
identical  between  Sij  and  s,-+Ufj+v.  If  the  two  summations  can  be  performed  simultaneously  (which 
is  possible,  using  the  architecture  described  in  Section  2.3),  then  accessing  a  single  additional 
composite  data  sample  permits  an  additional  integrator  output  to  be  computed.  This  idea  can  be 
carried  forward  so  that  each  additional  data  access  permits  an  additional  integrator  output  to  be 
computed. 

This  savings  can  be  generalized.  For  simplicity,  assume  that  u  >  0  and  v  >  0.  (The  algorithm 
can  be  applied  to  negative  velocities  as  well.)  Note  that 


TV+n-l 

*i+nu,j+nv  =  ^  Bk-n(ai+uk,j+vk)  •  (27) 

k=n 

For  each  i,  j  pair  a  set  of  integrator  outputs  can  be  computed  for  some  range  of  n.  One  needs  to 
consider  the  set  of  z,  j  indices  and  the  range  of  n.  The  set  of  z,j  is  easily  determined.  Assuming 
that  all  possible  values  of  n  are  to  be  used,  it  is  only  necessary  to  search  set  5, 

S  =  {(i,j);  0  <  i  <  u  or  0  <  j  <  v}  .  (28) 

This  set  is  illustrated  by  the  dotted  region  in  Figure  21.  Initial  target  positions  outside  this  set  are 
obtained  for  values  of  n  /  0.2  For  example,  the  black  circle  within  the  dotted  region  indicates  a 
particular  i,j  pair  with  n  =  0.  Values  of  n  >  0  correspond  to  initial  target  positions  indicated  by 
the  sequence  of  black  circles  in  the  clear  region. 

The  range  of  n  can  now  be  computed.  If  the  desire  is  to  detect  only  those  targets  that  do  not 
enter  or  exit  the  field  of  view  during  the  N  frames,  then  the  maximum  value  of  n  is  determined  by 
restricting  n  so  that  the  indices  i  +  uk,j  +  vk  in  Equation  (27)  satisfy 

0  <  i  +  uk  <  Mx  —  0.5 

0  <  j  +  vk  <  My  -  0.5  .  (29) 


In  this  case 


2This  statement  is  not  completely  accurate;  because  u  and  v  are  floating  point  numbers  and  nearest 
neighbor  interpolation  is  used,  it  is  possible  to  miss  some  initial  positions.  These  omissions  can  be 
eliminated  by  using  the  set  S  =  {(i,j);  0  <  i  <  u  or  0  <  j  <  v)  at  the  expense  of  duplicating  some 
initial  positions.  This  additional  expense  is  small  relative  to  the  overall  cost. 
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My-  1 
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Mx-  1 


0  <  n  <  min 


Figure  21.  Search  region  foriyj. 


(30) 


for  u  /  0  and  v  ^  0.  [If  either  u  =  0  or  v  =  0,  then  the  appropriate  fraction  is  deleted  from 
Equation  (30).]  The  number  of  integrator  outputs  computed  for  a  given  i,j  pair  is  then3 


r  .  /mx-o 
rni — u 


—  0.5  —  i  My  —  0.5  —  j 


)]  -  (JV  -  1)  . 


(31) 


The  number  of  data  accesses  [i.e.,  values  of  k  determined  from  the  bounds  in  Equation  (29)]  is 

r  .  / Mx  —  0.5  —  »  My- 0.5 -j\1  ,  x 

|mm( — : — .  )|  ’  <32> 


3The  notation  [x]  is  used  to  indicate  the  ceiling  of  x,  the  smallest  integer  that  is  greater  than  or 
equal  to  x.  Similarly,  [ x J  is  the  floor  of  x,  the  largest  integer  that  is  less  than  or  equal  to  x. 
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and  the  ratio  of  outputs  computed  per  data  access  is 


1  - 


_ TV  -  1 _ 

[min((Mx  -  0.5  -  i)/u,(My  -  0.5  -  j)/v) ] 


(33) 


This  ratio  approaches  unity  as  the  frame  size  becomes  much  larger  than  the  distance  a  target  moves 
from  one  frame  to  the  next. 

Alternatively,  the  range  of  n  may  be  chosen  to  include  those  targets  that  enter  or  exit  the  field 
of  view  during  the  TV  frames.  Admittedly,  such  targets  must  have  a  higher  SNR  to  be  detected.  In 
this  case  the  range  of  n  is  determined  by  the  bounds 


o  < 

i  +  u(N  +  n  —  1) 

(34) 

o  < 

j  +  v(N  +  n  —  1) 

(35) 

i  +  un 
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o 

bi 

(36) 

j  +  vn 

<  My  -  0.5  , 

(37) 

which  give  the  following  range  on  n. 


max 


(l  -  N  -  1  -  N  -  <  n  <  min  ( 


Mx  —  0.5  —  i  My  —  0.5  —  j' 


(38) 


Given  the  bounds  on  i,  j  and  noting  that  n  is  an  integer,  the  lower  bound  on  n  can  easily  be  shown 
to  be  1  —  TV  so  that  the  range  of  n  becomes 


1  —  TV  <  n  <  min 


( 


Mx  -  0.5  -  i 
u 


My  -0.5-j 
v 


) 


(39) 


The  number  of  integrator  outputs  computed  for  a  given  i,j  pair  is  the  number  of  values  of  n  in  the 
above  equation, 


|^min  f 


Mx  -  0.5  -  i  My  —  0.5  —  j' 


u 


+  N-  1 


(40) 


The  number  of  data  accesses  is  identical  to  the  previous  case.  The  ratio  of  outputs  per  data  access 
in  this  case  is  then 


1  + 


_ N  -  1 _ 

[min((Mx  -  0.5  -  i)/u,(My  -  0.5  -  j)/v) ] 


(41) 
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which  is  clearly  greater  than  unity. 

The  resulting  algorithm,  which  tests  a  given  velocity  hypothesis,  is  described  in  this  series. 

1.  Let  N  be  the  number  of  image  frames,  each  with  size  MXy  My  in  the  x,  y  dimensions. 

2.  Let  M  be  the  detection  threshold  (i.e.,  the  number  of  binary  quantized  frames  in 
which  a  target  must  appear  in  order  to  have  a  detection). 

3.  Let  a,  j  be  the  composite  of  the  binary  data  for  pixel  i,j  and  all  frames,  as  described 
in  Equation  (21). 

4.  Let  u,v  be  the  frame-to-frame  target  movement  in  x,j/  for  the  desired  velocity  hy¬ 
pothesis,  u  >  0  and  v  >  0. 

5.  For  each  pair  i,j  such  that  0  <  i  <  u  or  0  <  j  <  v,  do  the  following: 

(a)  Set  z*  =  0  for  k  =  1, ...  ,7V  —  1. 

(b)  Let  L  =  [min((Mx  —  0.5  —  i)/ti,  ( My  —  0.5  —  j)/v) ]  —  1. 

(c)  For  n  =  0, . . . ,  L,  do  the  following: 

i*  ^t+(n— N+i)ti,j+(n— N+i)v  =  l  "h  &N— l (^i+nuj+nv)*  Compare  this  result 

with  threshold  M . 

ii.  For  77T  —  A  1,...,2  set  Zm  —  Zm—\  -b  -Sm-l +nv ) • 
iii •  Set  Z\  =  50(at+nuj+nv). 

(d)  Forn  =  L  + 1, ...,  X  +  Ar  - 1,  the  outcome  is  5t-4.(n_jV+1)u>i+(n-/sr+i)t;  =  zN+L-n* 
Compare  this  result  with  threshold  M . 

2.3  Architecture  for  Binary  Integration 

An  architecture  that  implements  the  binary  integration  algorithm  is  shown  in  Figure  22.  For 
a  given  i,  j  pair,  this  algorithm  computes  the  st+nUfJ+nt,.  The  top  row  of  adders  and  registers  is  used 
to  implement  the  sums  and  partial  sums  zi, . . . ,  zjv-i  •  The  subtracter  compares  the  result  z/v-i 
with  the  threshold  A/.  Threshold  crossings  are  stored  in  a  FIFO  with  an  index  (counter  output). 

The  bottom  portion  of  Figure  22  generates  the  necessary  addresses.  The  register  stores  a 
base  address  corresponding  to  the  particular  i,j  pair  in  use.  The  counter  generates  the  sequence 
n  =  0,. . .  to  be  used  as  an  address  to  an  offsets  table,  which  is  loaded  from  an  external  processor 
and  contains  the  memory  offsets  corresponding  to  the  pixel  offset  nu,ni\  The  table  only  needs  to 
be  changed  when  the  velocity  vector  is  changed.  Finally,  the  adder  combines  the  base  address  with 
the  necessary  offset,  producing  the  memory  address. 

Operation  of  this  architecture  is  straightforward.  Initially,  the  top  row  of  registers  and  the 
FIFO  are  set  to  zero.  With  each  computation  cycle  the  counter  drives  the  address  generation 
portion  to  access  a  particular  word  of  memory.  The  composite  binary  sample  at  this  address  is 
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Figure  22.  Architecture  for  integration  step  of  algorithm. 


presented  to  the  inputs  These  bits  are  added  to  the  partial  sums  that  are  available 

from  the  register  outputs.  The  output  of  the  right-most  register  contains  the  completed  sum,  which 
is  compared  with  the  threshold  M .  At  the  end  of  each  computation  cycle,  the  adder  outputs  are 
stored  in  the  registers,  detections  are  stored  in  the  FIFO,  and  the  counter  is  incremented.  Step  5(d) 
of  the  algorithm  is  implemented  by  reading  N  —  1  words  of  zeros  from  memory.  These  zeros  can 
be  accessed  by  storing  appropriate  values  in  the  offset  table.4  Note  that  the  final  register  contents 
are  all  zeros  so  that  the  adder  pipeline  is  initialized  for  the  next  t,jf  pair. 

For  typical  parameter  values,  this  architecture  can  be  implemented  in  a  small  ASIC.  For 
example,  suppose  N  <  16,  MxMy  <  106,  and  max(Afr,  My)  <  4,000.  In  this  case  the  typical 
requirements  are  listed  in  Table  1.  A  24-bit  address  bus  is  assumed  for  the  memory.  The  FIFO 
depth  is  arbitrarily  set  to  256  words,  although  64  might  be  sufficient.  The  left-most  registers  and 
adders  require  fewer  bits  than  the  right-most.  It  is  clear  that  the  circuit  requirements  are  small. 

Note  that  the  architecture  in  Figure  22  is  easily  modified  to  permit  the  cascading  of  several 
such  ASICs  for  large  numbers  of  frames.  In  this  case  the  input  to  the  left-most  register  is  now  the 
output  of  an  adder  that  has  Bo  and  an  off-chip  datum  as  inputs.  In  addition,  the  output  of  the 
right-most  register  is  sent  off-chip.  Thus  the  top  portion  of  additional  ASICs  may  be  cascaded  with 
the  final  ASIC,  which  also  performs  detection  and  address  generation.  With  this  extension  scheme, 
the  limit  on  the  total  number  of  frames  is  set  by  the  width  of  the  registers,  adders,  and  FIFO. 

The  velocity  filter  architecture  is  shown  fitted  into  a  signal  processor  architecture  in  Fig¬ 
ure  23.  There  is  a  dual  port  memory  for  storing  the  data  prior  to  registration,  normalization,  and 


*High  speed  operation  may  require  using  pipeline  registers  at  the  address  output  and  the  Bk  inputs. 
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TABLE  1 

Architecture  Component  Requirements 


Item 

Amount 

Register 

Adder  (1  bit  -f  k  bit) 
Adder/Subtracter 

Counter 

RAM  (offset  table) 

FIFO 

78  bits 

49  bits 

29  bits 

12  bits 

4K  words  x  20  bits 

256  words  x  18  bits 

quantization.  A  smaller  memory  (a  single  frame’s  worth)  is  used  to  separately  store  the  binary 
quantized  data.  The  CPU  and  a  separate  memory  share  a  second  bus,  which  allows  the  CPU  to 
perform  other  computations  while  the  velocity  filter  is  in  operation. 


190269-20 


Figure  23.  Architecture  for  binary  integration. 


To  achieve  sufficient  processing  performance,  it  may  be  necessary  to  have  several  velocity 
filters  operating  in  parallel.  One  method  for  four  filters  (suggested  by  A.  H.  Huntoon)  is  shown  in 
Figure  24.  The  composite  memory  is  divided  into  a  number  of  submemories  equal  to  the  number 
of  velocity  filters.  Each  submemory  is  assigned  a  contiguous  set  of  columns  of  image  data.  The 
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Figure  24 •  Parallel  architecture  for  binary  integration. 


memories  are  connected  to  the  velocity  filters  through  a  data  exchange  circuit.  The  sequence  of 
tested  target  positions  can  be  arranged  so  that  during  each  memory  access  cycle,  each  velocity  filter 
needs  data  from  a  separate  submemory.  In  this  case  the  (i,j)  search  region  will  be  different  from 
that  depicted  in  Figure  21. 

2.4  Extending  the  Algorithm  to  Nonscanning  Sensors 

With  a  scanned  sensor,  a  moving  target  is  expected  to  move  less  than  one  pixel  during  the 
time  the  target  is  scanned  for  a  single  frame.  Thus  unresolved  targets  only  illuminate  a  single  pixel 
per  frame. 

Nonscanned  sensors  are  somewhat  different  because  all  pixels  are  collecting  data  during  the 
entire  frame  time.  The  result  is  that  in  addition  to  the  frame-to-frame  motion  expected  for  a 
moving  target,  the  target  moves  within  a  frame,  producing  a  streak  or  straight  line  within  each 
frame.  Thus  a  binary  integration  algorithm  must  sum  over  several  pixels  within  each  frame. 


33 


One  approach  to  summing  multiple  pixels  per  frame  is  to  implement  the  algorithm  and  ar¬ 
chitecture  described  above  to  account  for  the  frame-to-frame  motion,  then  combine  outputs  from 
the  velocity  filter  to  sum  over  multiple  pixels  per  streak.  This  combining  could  be  implemented  in 
an  ASIC,  but  is  probably  easier  to  implement  in  the  CPU  shown  in  Figure  23. 

This  algorithm  is  described  by  the  following  sequence.  The  basic  idea  is  to  form  a  one¬ 
dimensional  output  from  the  frame-to-frame  portion  of  the  velocity  filter,  which  can  then  be  con¬ 
volved  with  a  filter  representing  the  moving  target  within  a  frame.  The  algorithm  and  architecture 
described  earlier  are  used  to  obtain  the  one-dimensional  data  needed  for  the  convolution. 

1.  Let  u,v  be  the  frame-to-frame  movement  in  x,y  for  the  desired  velocity  hypothesis, 
u  >  0  and  v  >  0. 

2.  Let  6  =  arctan(n/u)  and  h  =  \/{u 2  +  v2).  The  angle  of  the  velocity  vector  is  0,  and 
h  is  the  length  of  that  vector.  The  terms  cos(0)  and  sin(0)  may  be  computed  here 
for  later  use. 

3.  Choose  an  initial  point  a,/3  for  the  integration,  similar  to  i,j  in  the  algorithm  de¬ 
scribed  in  Section  2.3.  In  this  case  it  is  necessary  to  have  either  a  =  0  or  (3  =  0,  and 
0  <  a  <  u  and  0  <  /3  <  v. 

4.  For  /  =  0, . . . ,  h  —  1,  do  the  following: 

(a)  Let  z  =  a  +  /  cos(0)  and  j  =  f3  +  /sin(0). 

(b)  Let  the  binary  detection  threshold  be  zero,  M  =  0. 

(c)  Implement  the  algorithm  described  in  Section  2.3,  using  the  architecture  in 
Figure  22  with  the  values  of  i,  j  and  M  specified  above. 

(d)  Outputs  Si+nuj+nv  are  saved  according  to:  f\+hn  =  $,-+nuIj+nv* 

5.  The  one-dimensional  data  set  fn  is  an  ordered  sequence  of  outputs  from  convolving 
the  frame-to-frame  portion  of  the  velocity  filter  with  the  data.  The  complete  velocity 
filter  output  is  computed  by  the  convolution  fn  *  pn,  where  gn  =  {1,...,1}  is  a 
sequence  of  h  ones  representing  the  target  streak  within  a  frame. 

6.  The  output  of  the  convolution  is  compared  with  the  desired  M-out-of-N  detection 
threshold. 

This  algorithm  may  be  viewed  as  postprocessing  the  result  of  the  algorithm  described  in  Sec¬ 
tions  2.2  and  2.3.  Postprocessing  is  a  convolution  by  a  rectangular  window  that  can  be  implemented 
efficiently  by  sliding  the  window  (adding  one  data  point  and  subtracting  another)  without  the  need 
for  an  FFT  approach. 
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3.  CONCLUSION 


A  maximum  likelihood  moving  target  detection  algorithm  is  derived  based  on  a  Gaussian 
model.  Data  samples  are  assumed  to  be  temporally  stationary,  implying  that  image  frames  are  reg¬ 
istered  and  that  there  is  no  moving  clutter.  Data  samples  are  allowed  to  be  spatially  nonstationary, 
corresponding  to  nonuniform  background  clutter. 

In  Section  1,  the  probability  distribution  of  normalized  experimental  data  is  shown  to  agree 
closely  with  theoretical  predictions.  This  agreement  holds  for  a  variety  of  clutter  backgrounds 
(provided  the  temporal  stationarity  requirement  is  met)  and  validates  studying  the  theoretical 
model  as  a  means  of  predicting  system  performance,  which  is  evaluated  in  terms  of  the  probabilities 
of  detection  and  false  alarm. 

In  some  system  applications  the  computation  and  memory  requirements  of  this  algorithm 
may  be  large,  particularly  when  either  real-time  processing  or  reduced  processor  size  is  required. 
A  binary  integration  algorithm  is  developed  in  Section  2  to  address  these  issues.  Analyzing  its 
detection  and  false  alarm  performance  shows  a  2-dB  loss  relative  to  the  full-precision  algorithm  that 
is  discussed  in  Section  1.  Similar  results  are  obtained  in  radar  applications.  To  efficiently  implement 
the  binary  integration  algorithm  an  architecture  concept  is  described  that  takes  advantage  of  a  small 
ASIC  to  perform  the  bulk  of  the  computation. 

One  of  many  areas  open  for  further  research  is  frame  registration;  in  many  system  appli¬ 
cations  the  sensor  is  moving,  requiring  registration  for  sequential  frames.  Potential  effects  to  be 
compensated  include  translation,  scaling,  and  geometric  distortion,  and  strongly  depend  on  sensor 
design.  Other  registration  issues  arise  due  to  moving  (or  apparently  moving)  clutter.  In  this  case 
the  entire  frame  is  not  adjusted;  instead  the  clutter  region  must  be  selectively  corrected.  Difficulties 
arise  along  the  edges  of  the  moving  clutter. 

Another  research  area  is  incorporating  target  statistics  models  into  the  detection  algorithm. 
The  maximum  likelihood  algorithm  described  in  Section  1  permits  a  single  data  sample  of  suffi¬ 
cient  amplitude  to  cause  a  detection.  An  algorithm  that  requires  a  more  distributed  target  signal 
component  might  offer  improved  performance. 

Adapting  the  maximum  likelihood  algorithm  to  time- varying  clutter  backgrounds  also  requires 
further  work.  Not  all  backgrounds  (e.g.,  sea  glint)  satisfy  the  temporal  stationarity  requirement. 
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APPENDIX  A 

MAXIMUM  LIKELIHOOD  ALGORITHM  DERIVATION 


A  maximum  likelihood  algorithm  (presented  briefly  in  an  earlier  paper  [3])  is  developed  for 
moving  target  detection.  The  following  provides  a  more  detailed  derivation,  using  the  model  that  (in 
the  absence  of  targets)  pixel  samples  are  independent  normal  random  variables  that  are  temporally 
stationary.  In  other  words,  the  underlying  mean  and  variance  of  any  given  pixel  is  assumed  constant 
from  frame  to  frame.  These  means  and  variances  are  considered  to  be  unknown  and  to  vary  spatially 
(i.e.,  from  pixel  to  pixel). 

This  model  is  based  on  independent  arrival  of  photons  in  each  detector  (a  Poisson  random 
variable  approximated  a s  normal).  The  average  photon  count  in  a  given  detector  corresponds  to 
the  image  intensity.  The  variance  in  the  photon  count  is  referred  to  as  “photon  noise.”  Additional 
sources  are  thermal  and  readout  noise,  which  are  also  assumed  to  be  independent  from  sample  to 
sample. 

Stationarity  over  time  implies  that  all  frames  in  a  sequence  are  spatially  registered  and  that 
there  is  no  moving  clutter.  Nonstationary  spatial  statistics  arise  due  to  the  nonuniform  intensity 
of  the  background  scene. 

Prior  to  the  algorithm  derivation,  it  is  useful  to  establish  a  notation  that  simplifies  the  equa¬ 
tions.  Unprocessed  image  data  are  represented  by  where  (i,j)  are  the  spatial  indices  and  ( k ) 
is  the  temporal  index.  The  set  of  indices  that  represent  a  given  target  motion  hypothesis  is  denoted 
by  S.  It  is  assumed  that  the  target  moves  a  minimum  of  one  pixel  between  frames.  Set  S  contains 
M  members  with  each  member  corresponding  to  a  different  pixel,  (t,jf).5 

A  simplified  notation  is  used  for  the  N  frames  of  the  M  pixels  in  5.  Pixel  coordinates  are 
mapped  to  a  single  sequential  index.  Additionally,  the  temporal  index  is  permuted  so  that  the 
hypothesized  target  is  in  those  samples  that  have  a  temporal  index  of  zero.  To  define  this  notation, 
let  the  pixels  in  S  be  represented  by  the  triplet  for  0  <<  AT.  Data  samples  from  these 

M  pixels  and  N  frames  are  represented  by  pmin,  where 

Pm,n  =  rt  mjmj  (A.l) 

/  =  (n  +  km)modN  (A. 2) 

and  0  <  m  <  M ,  0  <  n  <  N . 


5 For  a  scanning  sensor,  an  unresolved  target  occupies  one  pixel  per  frame,  giving  M  =  N .  For 
a  staring  sensor,  the  target  may  move  during  the  observation  time  and  occupy  several  pixels  per 
frame,  giving  M  >  N . 
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For  the  sample  pm}U,  the  unknown  mean  and  variance  of  the  noise  are  pm  and  respectively. 
The  unknown  target  component  of  pm, o  is  sm. 

The  goal  in  this  derivation  is  to  develop  a  likelihood  function  to  be  tested  in  choosing  between 
two  hypotheses.  Let  H\  be  the  hypothesis  that  a  target  is  present,  and  let  H o  be  the  hypothesis 
that  no  target  is  present.  The  probability  density  of  receiving  data  given  hypothesis  Hn,  is 

fp({pi,k}\Hn)  .  (A. 3) 

The  likelihood  function,  A,  is  defined  as  the  ratio  of  this  density  for  the  two  hypotheses,  where  the 
densities  are  maximized  over  the  unknown  parameters  for  clutter  and  signal  [4]. 

max{/i(fg||J(}  f p({pi,k}\H\) 
max{/i.,<r;}  f p{{pi,k}\Ho) 

The  numerator  is  maximized  over  the  set  of  /xt  ,  crt-,  and  st  for  0  <  i  <  M .  Similarly,  the  denominator 
is  maximized  over  the  set  of  fi{  and  <7,  for  0  <  i  <  M . 

The  msiximization  of  the  denominator  is  derived  first.  Because  the  p+j.  are  independent 
normal  random  variables,  the  density  fp  is 


M-l  N—l 


fP({pi,k}\Ho)  =  n  II  ~?L— exp(-(p,,fc  - /z,)2/2a?) 
isso  k=o  y27rat2 


(A.5) 


Taking  the  partial  derivative  of  this  function  with  respect  to  the  parameter  /x,  gives 


d 


N- 1 


Px,k  ~  Pi 


dp/p  fp  ^  -2 


(A.6) 


k=0 


where,  for  convenience,  the  argument  to  the  function  fp  is  dropped.  When  this  partial  derivative 
is  set  to  zero,  it  is  seen  that  fp  is  maximized  when  px  =  a,*,  where 

l  N~l 

a,  =  jv^  p,’k  •  (A-7) 

k= 0 

Similarly,  substituting  a;  for  and  taking  the  partial  derivative  of  fp  with  respect  to  the 
parameter  cr2  gives 
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(A.8) 


J_f  _fNyl(pi.k-ai)2-^ 

W  ph 

*  A;=0  1 


Setting  Equation  (A.8)  to  zero  shows  that  fp  is  maximized  for  <7t2  =  c,,  where 


N-i 


'<  =  i  !>.* - “O1 


k=0 


(A.9) 


Finding  the  values  of  /it,  07,  and  st  that  maximize  the  numerator  of  Equation  (A. 4)  is  more 
complex.  In  this  case  the  density  fp  is  given  by 


M—\  N- 1 


fP({p,,k) \Hi)  =  II  II  /  =  exP (~(Pi,k  -  Pi  -  Si6(k)f/2a?) 

1=0  k=Q  yJ27T<jf 


(A. 10) 


where  S(k)  is  the  Dirac  delta  function,  defined  to  be  unity  for  k  =  0,  and  zero  elsewhere.  Taking 
the  partial  derivative  of  fp  with  respect  to  these  three  parameters  yields  the  following  equations: 


fN- 1 


Pi,k  -  Pi  Si 


-rk  )h 


(A. 11) 


d_ 

daV 


Q  _2  fp  ~ 


( (p.,0  -  Pi  -  -s.)2  ~ a,2  iRiA  -ptf  - 


—f  = 

dsiJp 


V 

Pt,0  Ab'  5 


+  £ 

k=  1 


2a? 


/p 


(A.12) 


(A. 13) 


Setting  Equations  (A. 11),  (A.12),  and  (A, 13)  to  zero  and  solving  the  results  simultaneously  shows 
that  fp  is  maximized  for 


/it  =  Ui 


=  Vi 


Si  =  no¬ 


where 


Ui 


1 

N  -  1 


N- 1 


fc=l 


(A. 14) 
(A. 15) 
(A.16) 


(A. 17) 
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V.  =  4  £  (pi.*  -  «i)2  (A-18) 

1  k=  1 

Wi  =  p,)0  -  •  (A. 19) 

Note  that  Equations  (A. 17)  and  (A. 18)  are  different,  particularly  in  the  limits  of  the  summations, 
from  (A. 7)  and  (A. 9). 

Substituting  the  density  functions  in  Equations  (A. 5)  and  (A.  10)  into  Equation  (A. 4)  gives 
the  following  expression  for  the  likelihood  function  when  the  numerator  and  denominator  are  max¬ 
imized: 


M-i 


A  = 


mo 


n^o1 


v  N/2  Yl^'expi-ipi'k -u,y/2v) 


c~N/2  Uk=o  exP (~(pi,k  ~  a.)2/2c) 


It  is  often  easier  to  work  with  the  logarithm  of  the  likelihood  function,  which  is 
£  =  ln(A) 

N-l  (  \2  AT— 1 

( Pi*  ~ 

k=l  2v{  fc=0 


M— 1 

=  £ 

t=0  L 


-N 


In 


^  _  'g  (Pi,*-**,)2  +  £  (Pi,*  -  a.)2 


2  \CiJ  f-t  2vi  '  2c, 

This  equation  is  simplified  by  first  noting  from  Equations  (A. 9)  and  (A. 18)  that 


N-i 


(p«,*  ~  Q«)2  _  {Ehb.  —  u')2  —  o 


*=o 


*=i 


2vx 


Further  simplification  is  obtained  using  the  approximation 
ln(l  +  x)  ^  x 

for  |x|  <C  1-  The  logarithm  term  in  Equation  (A.21)  is  then  approximated  by 

1  (pt,o  ~~  Qt)2  \ 


In 


(5)- 


In  1  - 


N-i  a 


-1  (p,,o-Qj)2 

N-i  a 


(A. 20) 


(A.21) 
(A. 22) 


(A. 23) 


(A. 24) 


(A. 25) 

(A.26) 
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In  evaluating  the  accuracy  of  this  approximation,  note  that  s ?  =  (pti o  —  a**)2  is  an  estimate  of  the 
target  power  in  pixel  i,  and  a?  =  c;  is  an  estimate  of  the  noise  power.  Thus  the  approximation 
is  valid  when  the  SNR  is  much  smaller  than  the  number  of  frames  processed.6  The  resulting 
log-likelihood  function  is 


£  = 


N  (/>;, o  -  a,)2 

-  U  ^  ~ 


2 (N  -  1)  « 

The  resulting  maximum  likelihood  detection  algorithm  can  now  be  stated  as 


(A.27) 


fccj 

*-* 

IV 

t-fc. 

(A.28) 

Hq  \i  C  <  t  , 

(A. 29) 

where  t  is  the  detection  threshold  (often  set  for  a  constant  false  alarm  rate  criterion)  and  where, 
using  the  original  notation  for  unprocessed  data  samples, 

C-  =  £ 

(i,j,k)eS 

(A. 30) 

J  _  —  Pi.j)2 

-  d2  . 

(A.31) 

l  N-l 

fa  =  jf  E  ru* 

k=0 

(A. 32) 

1 

k=  0 

(A. 33) 

Set  5  is  the  set  of  spatial-temporal  indices,  (t,j,fc),  corresponding  to  the  target  motion  hypothesis 
H\.  The  constant  multiplier  in  Equation  (A.27)  is  dropped,  noting  that  the  constant  can  be 
incorporated  into  the  detection  threshold.  The  term  ft is  the  sample  mean  of  N  temporal  samples 
of  pixel  (i,  j).  The  term  afj  is  the  sample  variance  of  the  same  N  temporal  samples  of  pixel  (i,j). 


6In  the  case  where  the  SNR  is  not  much  smaller  than  the  number  of  frames,  the  approximation  to 
the  log-likelihood  function  is  less  accurate.  In  such  cases,  however,  the  SNR  is  large  and  some  loss 
can  be  tolerated. 
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The  term  is  referred  to  as  the  “normalized”  data.  Finally,  for  this  detection  algorithm  to  be 
meaningful  the  number  of  frames  processed  must  satisfy  N  >  3. 
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APPENDIX  B 

SINGLE  PIXEL  PROBABILITY  DISTRIBUTIONS 


The  normalized  data  in  the  maximum  likelihood  electro-optic  moving  target  detection  algo¬ 
rithm  is  shown  to  have  a  beta  distribution,  which  is  used  to  determine  the  detection  threshold, 
probability  of  false  alarm,  and  probability  of  detection  for  both  the  full  precision  maximum  likeli¬ 
hood  detection  and  the  binary  integrator. 

The  log-likelihood  function  C  of  the  maximum  likelihood  algorithm  is  expressed  a s  the  sum 
of  normalized  data  where 


c  = 

X]  di,j,k 

(' iJ.*)€S 

(B.l) 

di,j,k  = 

*.2,; 

,  N- 1 

(B.2) 

c*. 

II 

n= 0 

j  N-l 

(B.3) 

dh  = 

Tr  —  W»i)2  ’ 

n= 0 

(B-4) 

and  5  is  the  set  of  space  (i,j)  and  time  ( k )  indices  corresponding  to  the  target  motion  hypothesis. 

The  probability  distribution  fk(x|iV)  of  is  shown  to  follow  a  beta  distribution  for  the 
noise-only  case  and  a  noncentral  beta  distribution  when  a  target  is  present.  For  the  radar  case, 
where  fi  is  known  a  priori  to  be  zero  and  where  is  complex,  this  distribution  is  well  known  [5]. 

B.l  Distribution  for  Noise  Only  (No  Target) 

The  direction  taken  is  to  relate  the  distribution  i^(x|A)  of  to  the  beta  distribution.  For 
notational  simplicity,  the  pixel  indices  (i,j)  will  be  dropped  because  they  never  change  during  the 
derivation. 

The  beta  probability  density  function  (pdf)  with  parameters  a  and  / 3  is 


fB(x\a,f3) 


T(a  +  /?)  o—i 

T(a)TU3) 


(1  -  x)«-> 


(B.5) 


for  0  <  x  <  1,  and  /s(x|a,/3)  =  0  elsewhere.  If  Xi  and  X2  are  independent  chi-square  random 
variables  with  v\  and  v 2  degrees  of  freedom  respectively,  then  the  ratio 
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(B.6) 


Xi 

Xi  +  x2 

has  the  beta  distribution 


Fb 


2  ’  2  ) 


(B.7) 


The  random  variable  djt  can  be  written  in  a  form  similar  to  Equation  (B.6)  by  using  a  linear 
transformation  of  independent  random  variables.  Let  Q  be  an  orthogonal  TV  x  TV  matrix,  and  let  & 
be  a  length-TV  vector  of  independent  standard  normal  random  variables  (i.e.,  zero  mean  and  unity 
variance),  where  the  nth  element  is 


rn  -  M 

Un  =  - 

a 

for  0  <  n  <  TV.  Then  the  vector 
v  =  Qu 


(B.8) 


(B.9) 


is  also  a  vector  of  independent  standard 


normal  random  variables  vn,  where  0  <  n  <  TV. 


Many  possible  orthogonal  matrices 


can  be  constructed.  Let  the  first  and  second  rows  of  Q  be 


(B.10) 


and 


N 

N  -  1 


(B.ll) 


where  e*.  is  a  vector  with  a  one  in  the  Arth  position  and  zeros  elsewhere.  It  is  easy  to  verify  that 
these  two  rows  are  orthogonal,  and  each  has  unity  norm  as  is  necessary  for  Q  to  be  orthogonal. 

Let  u  represent  the  sample  average  of  the  unj 


1  7V”1 

« =  jj  E •  (B-12) 
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It  is  easy  to  show  that 


Vo 

=  uy/N 

(B.13) 

Vl 

-  1 

(B.14) 

TV-1 

TV-1 

Yrt 

n= 0 

=  Y<  . 

n=0 

(B.15) 

The  transformation  of  independent  standard  normal  random  variables  is  next  used  to  derive 
the  distribution  of  d Note  that  ( r —  ft)  is  related  to  un  by 


(rn  -  /*)  =  (tin  -  u)< 7 


(B.16) 


This  equality,  along  with  Equations  (B.13)  and  (B.15),  is  used  to  give  the  following  relation  between 
the  denominator  in  the  definition  of  d in  Equation  (B.2)  and  the  vn: 


TV- 1 


E(r"  ~  £)2 


TV— 1 

Y  ”  “)2 

n= 0 


TV-1 


=  CT2  - Nu 2  +  Y  Ur 


n= 0 


TV-1 


=  <72 


Y< 


(B.17) 


n=l 


Combining  Equations  (B.2),  (B.4),  (B.14),  and  (B.17)  gives  the  following  representation  for  df. : 

(Tk  -  A)2 


dk  = 


N 


»? +z^‘^<JV  x) 


(B.18) 

(B.19) 


Because  the  vn  are  independent  and  have  standard  normal  distribution,  the  term  u2  is  chi-square 
with  one  degree  of  freedom,  and  J2n=2  vn  is  chi-square  with  N  —  2  degrees  of  freedom.  Referring 
to  Equation  (B.6),  it  is  clear  that  dk  has  the  beta  distribution 
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Fd(x\N)  =  Fb 


x 


(B.20) 


N  -  1 


1  N  —  2 
2’  2 


where 


0  <  x  <  N  —  1 


(B.21) 


B.2  Distribution  with  Target 

The  probability  distribution  of  a  normalized  pixel  with  a  target  present  is  shown  to  be  a 
noncentral  beta  distribution.  In  this  case  the  chi-square  random  variable  in  the  numerator  of 
Equation  (B.6)  is  derived  from  normal  distributed  random  variables  with  nonzero  means.  As  in 
Section  B.l,  a  linear  transformation  is  applied  to  the  data  to  obtain  the  desired  form. 

Let  s  be  the  target  component  of  the  data  sample  r and  define  n  to  be  a  length  N  vector 
formed  from  the  data  samples,  where 


un  - 


Tn-  H 


Vs 


(B.22) 


These  un  are  independent  norma]  random  variables  with  unity  variance  and  nonzero  mean. 
Let  Q  be  the  orthogonal  matrix  defined  in  Section  B.l.  Then  the  vector 


21  =  Qu 


(B.23) 


is  also  a  vector  of  independent  normal  random  variables  with  unity  variance  and  nonzero  mean. 
The  derivation  outlined  in  Equations  (B.12)  through  (B.19)  is  repeated  to  show  again  that 


dk  = 


v\  +  £*=2  vl 


<N-  1) 


(B.24) 


It  will  be  shown  that  Vi  is  a  noncentra]  chi-square  random  variable  with  one  degree  of  freedom,  and 
Yln=2  vn  1S  a  central  chi-square  random  variable  with  N  —  2  degrees  of  freedom.  The  normalized 
data  dk  then  have  a  noncentra]  beta  distribution. 


From  Equations  (B.13)  and  (B.14),  the  means  of  the  random  variables  uo  and  uj  are 
E{v  o}  =  0 


(B.25) 
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E{vi } 


s 


(B.26) 


a 


IN-1 


N 


The  means  of  the  remaining  un  are  shown  to  be  zero  by  evaluating  the  sum  of  the  squares  of  the 
means  of  the  un,  using  Equations  (B.25)  and  (B.26). 


N- 1 

E  £K}2 

AT-1 

= 

(B.27) 

n= 1 

n=0 

TV— 1 

(B.28) 

n=0 

s2  N  -  1 

(B.29) 

(7 2  N 

= 

(B.30) 

Subtracting  E{v i}2  from  both  sides  of  Equation  (B.30)  proves  that  only  V\  has  nonzero  mean. 
Thus  Equation  (B.24)  indicates  that  dk  has  the  noncentral  beta  distribution  given  by 


Fi(x|JV,A)  =  Fc 


1  N  —  2 
2’  2 


(B.31) 


where 


0  <  x  <  N  -  1 


(B.32) 


and  A  is  the  noncentrality  parameter 


A  = 


s2  N  -  1 
C72  N 


(B.33) 


The  noncentral  beta  distribution  can  be  written  in  terms  of  the  central  beta  distribution 
according  to 


Fb'(x\oi,P,\) 


£e-V2(V2Tf B(x\a  +  „,p) 

n  U- 

n= 0 


(B.34) 


This  formula  is  easily  derived  from  one  that  is  similar  for  the  noncentral  F  distribution  ([6],  p.  192). 
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It  is  interesting  to  note  that  for  a  given  number  of  image  frames  TV,  the  noncentrality  parame¬ 
ter  A  is  the  product  of  a  constant  with  the  signal-power  to  noise-power  ratio  s2/a 2,  which  is  referred 
to  as  an  “instantaneous”  signal-to-noise  ratio,  SNRj.  A  similar  dependency  occurs  for  many  other 
detection  problems.  When  the  signal  power  is  zero,  the  noncentral  distribution  in  Equation  (B.34) 
reduces  to  a  single  central  beta  distribution,  as  it  should. 
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APPENDIX  C 

LIKELIHOOD  FUNCTION  STATISTICS 


The  probability  distribution  of  the  log-likelihood  function  is  needed  to  determine  the  proba¬ 
bilities  of  detection  and  false  alarm  for  the  maximum  likelihood  algorithm  and  is  expressed  in  terms 
of  the  density  function  derived  in  Appendix  B  for  normalized  pixels. 

The  log-likelihood  function  probability  distribution  Fc(x)  depends  on  the  SNR  through  the 
parameter  A  =  ^^SNRj  in  the  single  pixel  distribution  function,  F(i(x|Ar,  A).  In  computing  the 
probability  of  false  alarm,  the  SNRj  is  zero.  In  computing  the  probability  of  detection,  the  statistical 
distribution  of  SNRj  should  be  taken  into  account  in  computing  Fc(x).  The  distribution  of  SNRj 
depends  on  the  statistics  of  receiving  target  photons  as  well  as  those  of  background  clutter  intensity. 
The  approach  taken  here  is  to  compute  A)  for  constant  SNRj.  The  results  are  useful  in 

discussing  target  detection  in  terms  of  an  effective  signal-to-noise  ratio,  SNRefp,  as  described  at  the 
end  of  Section  1.4. 

The  log-likelihood  function  C  is  the  sum  of  normalized  data  d{jt 

£=  £  diij<k  ,  (C.l) 

(itjyk)eS 


where  S  is  the  set  of  coordinates  for  the  target  velocity  hypothesis.  The  model  used  in  the  maximum 
likelihood  derivation  in  Appendix  A  implies  that  are  independent  random  variables;  therefore, 
the  density  function  fc{x)  of  C  is  the  convolution  of  the  densities  of  The  single  pixel  density 

function  is  obtained  by  differentiating  Equation  (B.31)  with  respect  to  x,  yielding 


fd(x\N,X) 


1 

N  -  1 


x 

N  -  1 


1 

2’ 


N  -  2 
2 


(C.2) 


where  fgi  is  the  noncentral  beta  density. 

Let  M  be  the  number  of  elements  in  the  set  S  (i.e.,  M  is  the  number  of  samples  in  the  target 
motion  hypothesis).  The  log-likelihood  density  function  is  then  obtained  as  the  convolution  of  M 
density  functions, 


fc(x\N ,X)  =  fd  *  •  •  •  *  fd(x\N,  A) 


(C.3) 


The  distribution  function  is  the  integral  of  the  density  function,  which  can  also  be  expressed  as  a 
convolution  of  density  functions  with  a  unit  step  function,  h(x), 


Fc(x\NtX)  =  f  fc(u\N,X) 

J  —  OO 


du 


(C.4) 
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h*fc(x\N,X) 

(C.5) 

h  *  fd  *  *  *  *  *  fd(x |iV,  A)  , 

(C.6) 

where 


h{x) 


0,  x  <  0 
1,  x  >  0 


(C.7) 


The  probabilities  of  detection  and  false  alarm  are  the  complement,  1  —  Fc(x\ Ar,  A),  of  the 
distribution  function  and  may  be  expressed  as  the  convolutions 


Pd(x\N,\)  =  (1  -  /i)  */<**•••  *  fd(x\N,  A) 


(C.8) 


=  (1  -  h)  *  fd  *  •  •  -  *  fd(x\N,  A  =  0) 


(C.9) 


where  x  is  the  detection  threshold.  Simple  formulas  for  these  probabilities  are  difficult  to  obtain. 
Convolutions  can,  however,  be  evaluated  through  numerical  integration;  the  remainder  of  this 
appendix  discusses  the  details.  Results  are  plotted  in  Figure  11  where  N  =  10-frame  processing  for 
a  scanning  sensor  (i.e.,  M  =  N ). 

The  Af-fold  convolution  in  Equation  (C.8)  can  be  computed  as  an  Mth  order  integral.  The 
difficulty  in  such  an  integration  is  maintaining  accuracy  while  computing  numerical  results  in 
a  reasonable  time.  The  approach  taken  for  each  convolution  is  to  sample  the  function  to  be 
integrated  at  uniform  intervals  along  the  x-axis  and  use  Simpson’s  rule  [7,  8]  to  compute  the 
resulting  convolution  at  these  same  sample  points.  This  approach  is  applied  M  times  in  succession 
to  obtain  the  final  convolution. 

The  first  convolution  in  Equation  (C.8)  is  easily  shown  to  be 

(1  -  h)  *  fd(x\N,\)  =  l-Fd(x\N,\)  (C.10) 


=  Fq> 


( i  !Lzl  I  \\ 

V  N-  1  2  ’  2’  / 


(C.ll) 


and  can  be  computed  without  using  numerical  integration.  This  function,  and  each  subsequent 
convolution  of  fd(x\N,  A)  with  this  function,  is  the  complement  of  a  probability  distribution. 

Let  K ,  a  positive  even  number,  be  the  number  of  sample  intervals  in  the  range  0  <  x  <  N  —  1, 
and  let  6  be  the  sample  interval, 


6  = 


N  -  1 
K 


(C.12) 
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For  N  >  3  the  convolutions  are  computed  by  first  defining  two  sample  sequences,  gn  and  Un\ 
where 


9n 


'  fdW\N,\),  n  =  0 
<  fd(n6\N,  A),  n  =  1 , . . . ,  A' 
0,  n  >  K 


(C.13) 


f 

1, 

<  Fb*  (l 


71  <  0 

2’^)’  n  =  !»•••»  K 
n>  K 


(C.14) 


The  superscript  (m)  on  indicates  the  number  of  convolutions  represented  by  the  sequence. 
Convolutions  are  computed  using  gn  and  u ^  to  produce  tin™"1"1). 


The  parameter  7  is  carefully  chosen  in  the  range  0  <  7  <  1  to  avoid  a  singularity  in  fd(x\N,  A). 
[For  N  =  3  the  density  fd(x |  iV,  A)  also  has  a  singularity  at  x  =  2.]  As  x  approaches  zero  from  above, 
this  function  approaches  infinity.  The  parameter  7  is  chosen  so  that  when  Simpson’s  rule  is  applied 
directly  to  the  sequence  <7n,  representing  an  approximation  to  integrating  the  density  function 
fd(z | AT,  A),  the  result  is  very  close  to  unity.  For  example,  with  N  =  10,  A  =  0,  K  =  1,000,  and 
7  =  0.0701267113512541,  the  result  of  applying  Simpson’s  rule  to  gn  is  1  +  e,  where  \e\  <  10~16. 

Higher  order  convolutions  are  computed  according  to  the  formula 


u(m+1) 

uk 


1,  k  <  0 

<  |  Wn-k+KU^dk-n,  fc  =  1, .  •  •  ,  (m  +  1  )K 

0,  k  >  (m  +  1)K 


(C.  15) 


for  m  =  1, . . . ,  M  —  1,  and  where  the  weighting  sequence  wn  used  in  Simpson’s  rule  is  given  by 


min(mA',  k)  —  (k  —  K ) 

Wo,...,  W,nin(roK’,Jfc)-(Jfc-A') 

0 

{1} 

1 

{4,1} 

even 

{1,4,2,4,... ,4,2,4, 1} 

odd 

{4,2,4,. ..,4,2,4,!} 
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Simpson’s  rule  requires  an  even  number  of  sampling  intervals  (or  an  odd  number  of  sampling 
points).  In  the  case  where  the  number  of  intervals,  min(mA', k)  —  (k  —  /if),  is  odd,  the  additional 
term 


(c.i7) 

can  be  added  to  the  summation  in  Equation  (C.15)  for  iu_i  =  1;  however,  gx+i  =  0  and 
0  <  4-Li  —  1  so  that  this  additional  term  is  known  to  be  zero  and  can  be  omitted  from 
the  summation.  Accuracy  of  the  result  when  use  of  this  additional  term  is  implied  is  ensured  by 
noting  that  represents  samples  of  a  continuous  function  for  all  n,  and  gn  represents  samples 
of  a  continuous  function  for  n  >  1  (when  N  >  3). 

The  probability  of  detection  is  approximated  by 

PD(n6\N,X)  *  uiM)  .  (C.18) 

The  probability  of  false  alarm  is  obtained  with  the  same  formula  for  the  case  A  =  0. 
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